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Abstract. We study the Landau level spectrum using a multi-band k • p theory in 
monolayer transition metal dichalcogenide semiconductors. We find that in a wide 
magnetic field range the Landau levels can be characterized by a harmonic oscillator 
spectrum and a linear-in-magnetic field term which describes the valley degeneracy 
breaking. The effect of the non-parabolicity of the band-dispersion on the Landau 
level spectrum is also discussed. Motivated by recent magnetotransport experiments, 
we use the self-consistent Born approximation and the Kubo formalism to calculate the 
Shubnikov de-Haas oscillations of the longitudinal conductivity. We investigate how 
the doping level, the spin-splitting of the bands and the broken valley degeneracy of 
the Landau levels affect the magnetoconductance oscillations. We consider monolayer 
M 0 S 2 and WSe 2 as concrete examples and compare the results of numerical calculations 
and an analytical formula which is valid in the semiclassical regime. Finally, we briefly 
analyze the recent experimental results (Reference m) using the theoretical approach 
we have developed. 


PAGS numbers: 

1. Introduction 

Atomically thin transition metal dichalcogenides semiconductors (TMDCs) [H El E] 
are recognized as a material system which, due to its finite band gap, may have a 
complementary functionality to graphene, the best known member of the family of 
atomically thin materials. The experimental evidence that TMDCs become direct band 
gap materials in the monolayer limit [1] and that the valley degree of freedom [5] can be 
directly addressed by optical means PEIEIE] have spurred a feverish research activity 
into the optical properties of these materials uni ni EH E3]. Equally influential has 
proved to be the fabrication of transistors based on monolayer M 0 S 2 [II] which motivated 


LLs and SdH oscillations in monolayer TMDCs 


2 


a lot of subsequent research to understand the transport properties of these materials. 
Achieving good Ohmic contact to monolayer TMDCs is still challenging and this 
complicates the investigation of intrinsic properties through transport measurements. 
Nevertheless, signihcant progress has been made recently in reducing the contact 
resistance by e.g., using local gating techniques na, phase engineering lisi. making use 
of monolayer graphene as electrical contact imiHmn], or selective etching procedure 

ffl. 

Our main interest here is to study magnetotransport properties of monolayer 
TMDCs. Unfortunately, the relatively strong disorder in monolayer TMDC samples 
have to-date hindered the observation of the quantum Hall effect. Nevertheless, the 
classical Hall conductance has been measured in a number of experiments [151 CHI IHl 
l2^ EH] and was used to determine the charge density Ue and to extract the Hall mobility 
Hh- In addition, three recent works have reported very promising progress in the efforts 
to uncover magnetic held induced quantum effects in monolayer TMDCs. Firstly, in 
Reference m the weak-localization effect was observed in monolayer M 0 S 2 . Secondly, 
it was shown that in boron-nitride encapsulated mono- and few layer M 0 S 2 [H] and in 
few layer WSe 2 izni it was possible to measure the Shubnikov-de Haas (SdH) oscillations 
of the longitudinal resistance. Both of these developments are very signihcant and 
can provide complementary informations: the weak localization corrections about the 
coherence length and spin relaxation processes [251 ES], whereas SdH oscillations about 
the cross-sectional area of the Fermi surface and the ehective mass of the carriers. 

Here we hrst briehy review the most important steps to calculate the LL spectrum 
in monolayer semiconductor TMDCs in perpendicular magnetic held using a multi-band 
k ■ p model [3|. We show that for magnetic helds of R < 20T a simple approximation 
can be applied to capture all the salient features of the LL spectrum. Motivated by 
recent experiments in M 0 S 2 [H] and WSe 2 [20], we use the LL spectrum and the 
self-consistent Born approximation (SCBA) to calculate the SdH oscillations of the 
longitudinal conductance We discuss how the intrinsic spin-orbit coupling and the 
valley degeneracy breaking (VDB) of the magnetic held ahect the magnetoconductance 
oscillations. We also point out the diherent scenarios that can occur depending on the 
doping level. 

2. Landau levels in monolayer TMDCs 

Electronic states in the K and —K valleys are related by time reversal symmetry in 
monolayer TMDCs and hence in the presence of a magnetic held their degeneracy should 
be lifted. (Note that in the case of graphene the inversion symmetry, which is present 
there but not in monolayer TMDCs, ensures that in the non-interacting limit the LLs 
remain degenerate in the K and —K valleys.) Recently several works have calculated the 
Landau level (LL) spectrum of monolayer TMDCs using the tight-binding (TB) method 
[271 EH ES] and found that the magnetic held can indeed lift the degeneracy of the LLs in 
diherent valleys. However, due to the relatively large number of atomic orbitals that is 
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needed to capture the zero magnetic field band structure, for certain problems, such as 
the SdH oscillations of longitudinal conductance, the TB methodology does not offer a 
convenient starting point. On the other hand, a simplified two-band k-p model was used 
to predict unconventional quantum Hall effect [30] and to discuss valley polarization [31] 
and magneto-optical properties [32]. This model, however, did not capture the valley 
degeneracy breaking and was therefore in contradiction with the TB results and the 
considerations based on symmetry arguments. 

We first show that the VDB in perpendicular magnetic field can be described by 
starting from a more general, seven-bands k ■ p model [3]. To this end we introduce an 
extended two-band continuum model which can be easily compared to previous works 
[30l [311 [32l [33]. We then show a relatively simple approximation for the LL energies 
which will prove to be useful for the calculation of the SdH oscillations in Section [3l 


2.1. LLs from an extended two-band model 


Our starting point to discuss the magnetic field effects in monolayer TMDCs is a seven- 
band k ■ p model (fourteen-band, if the spin degree is also taken into account), we 
refer the reader to Refs. [3] for details. In order to take into account the effects of 
a perpendicular magnetic field, one may use the Kohn-Luttinger prescription, i.e., we 
replace the wavenumbers q = {q^, Qy) appearing in the seven-band model with operators: 
q — )■ q = j V -|- |A, where = (0, B^x, 0) is the vector potential in Landau gauge 
and e > 0 is the magnitude of the electron charge. Note that due to this replacement 
= Qx + iqy and = Qx — iqy become non-commuting operators: [g_, g+] = where 
\Bz\ is the strength of the magnetic field and [...] denotes the commutator. Working 
with a seven-band model is not very convenient and therefore one may want to obtain an 
effective model that involves fewer bands. This can be done using Lowdin-partitioning 
to project out those degrees of freedom from the seven-band Hamiltonian that are far 
from the Fermi energy. Since and g_ are non-commuting operators, it is important to 
keep their order when one performs the Lowdin-partitioning. To illustrate this point we 
first consider a two-band model (four-band including spin) which involves the valence 
and the conduction bands (VB and CB). We will follow the notation used in Reference 
[3]. One finds that the low-energy effective Hamiltonian in a perpendicular magnetic 
field is given by 


rj'Tj.S _ TT I UT,S I rjTiS 

-^eff — Bo + -h Hj^.p 

where s = 1 (s = — 1) denotes spin j' (j,) and 

_ q+q_ + 1 

Bo „ „ + ^QcI^bBzSz 

Imp 2 2 


( 1 ) 


( 2 ) 


is the free electron term [gp. 
Furthermore, 


Hlf = 


rAvbS^ 0 
0 rAcbS 


2 is the g-factor and is the Bohr magneton). 

(3) 
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describes the spin-orbit coupling in VB and CB [sz is a spin Pauli matrix) and r 
for the ±iP valleys. The k ■ p Hamiltonian I'sS'ds 


trTj.s _ I _i_ I tr''"’'* 

“kp ~ -^D “T-f^as “T -^310 “T -^cubJ 


±1 

(4) 


where 


H 


D 


£vh T 7r,sgh 
£cb 


(5a) 


rjTjS _ 

“as 




ar,sqlq- 


0 


_ __ 

■^cw 6 ,l Q 


0 ; ’ 

0 Kr,s{qlf \ 

0 ) ’ 

0 


r 


2 V ivi]h*qlqlq- + {v't1iyq-qir+ 


(2)5 


r]i]lqlq-q- + 7T^ighghg+ 


(56) 
(5 c) 

(sV) 


Here the operator gj is dehned as gj = g^; ± irqy. The material specihc properties are 
encoded in the parameters Svb, ^cb (band-edge energies in the absence of SOC), 7 t-,s 
( coupling between the VB and the CB) and a,-,®, /dT,s, «t,s, rj^l-, which describe the 
effects of virtual transitions between the VB (CB) and the other bands in the seven- 
band model. In general, the off-diagonal material parameters 'ys,T, and are 

complex numbers such that for the —K valley (r = —1) they are the complex conjugates 
of the K valley case (r = 1). In the absence of a magnetic held, the material parameters 
appearing in Eqs. (|5^ - (1^ can be obtained by, e.g., htting the eigenvalues of 
to the band structure obtained from density functional theory (DFT) calculations. We 
refer to Reference [3] for the details of this htting procedure and for tables containing 
the extracted parameters for monolayer semiconductor TMDCs. Here we only mention 
that such a htting procedure yields real numbers which depend on the spin index s but 
do not depend explicitly on the valley index r. (The parameters and cannot 
be obtained separately from htting the DFT band structure, only their sum, can 
be extracted. Fortunately, as we will see below, the ehect of ^ is very small in the 
magnetic held range we are primarily interested in. ) 

We note that a k - p model, similar to ours, was recently used in References [2^ 133] 
to calculate the LL spectrum. There are two diherences between our k ■ p Hamiltonian 
Eqs. dl]) and the model in References [29l [33]. The hrst one is that higher order terms 
that would correspond to our and ^ were not considered in References [29l [33] . 
We keep these terms in order to see more clearly the magnetic held range where the 
approximation discussed in Section 12.21 is valid. The second diherence can be found 
in our Hlf (1^ and the corresponding Hamiltonian used in [29l [33] . This diherence 
can be traced back to the way the magnetic held is taken into account in the ehective 
models that are obtained from multi-band Hamiltonians. In References [29l [33] hrst an 
ehective zero held two-band model was derived and then in a second step the Luttinger- 
prescription was performed in this ehective model. Therefore the terms which are 
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~ in the zero field case become ~ Q+Q- + Q-Q+ after the Luttinger-prescription. 
In contrast, as mentioned above, we perform the Luttinger prescription in the multi¬ 
band Hamiltonian and obtain the effective two-band model ([1]) in the second step. 
The two approaches may lead to different results because the operators g+, do not 
commute and this should be taken into account in the Lowdin-partitioning which yields 
the effective two-band model. 

The spectrum of can be calculated numerically using harmonic oscillator 
eigenfunctions as basis states. Taking 5^ > 0 for concreteness, one can see that the 
operators a and al defined as g_ = ^a, where Ib = \/h/{e\Bz\), satisfy 

the bosonic sbmmutation relation [a, = 1. (For < 0 one has to define = ^a, 

q_ = Therefore one can calculate the matrix elements of in a large, but 

finite harmonic oscillator basis and diagonalize the resulting matrix. For a large enough 
number of basis states the lowest eigenvalues of will not depend on the exact 
number of the basis states. Such a LL calculation is shown in Figured] for M 0 S 2 and 
in Figure d] for WSe 2 (we have used the material parameters given in Reference [5]). 
One can see that the LLs in different valleys are not degenerate and that the magnitude 
of the valley degeneracy breaking is different in the VB and CB and for the lower and 
higher-in-energy spin-split bands. While the results in the VB are qualitatively similar 
for M 0 S 2 and WSe 2 , considering the CB, for M 0 S 2 the valley splitting of the LLs is 
smaller in the higher spin-split band, whereas the opposite is true for WSe 2 . This is a 
consequence of the interplay of the Zeeman term in Eq. (|2]) and other, band-structure 
related terms which lead to VDB. (For M 0 S 2 the valley splitting in the higher spin-split 
CB (purple and cyan lines) is very small for the material parameter set used in these 
calculations and can only be noticed for large magnetic fields.) One can also observe 
that in the CB the lowest LL is in valley K, whereas in the VB it is in valley —K. 



Figure 1. Numerically calculated LL spectrum of M 0 S 2 . a) The first few LL in 
the higher spin-split VB. Red lines: the K valley (r = 1), blue lines: the —K valley 
(r = —1). The inset shows the LLs in the lower spin-split VB. b) The first few LL in 
the CB. LLs both in lower spin-split band and in the higher spin-split band are shown. 
Red and purple lines: the K valley, blue and cyan:the —K valley. 
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Figure 2. Numerically calculated LL spectrum of WSe 2 . a) The first few LL in 
the higher spin-split VB. Red lines: the K valley (r = 1), blue lines: the —K valley 
(r = —1). The inset shows the LLs in the lower spin-split VB. b) The first few LL in 
the CB. LLs both in lower spin-split band and in the higher spin-split band are shown. 
Red and purple lines: the K valley, blue and cyan:the —K valley. 


Further details of the VDB, including its dependence on the parameter set that can 
be extracted from DFT calculations, will be discussed in Section 12.21 Here we point out 
that these results qualitatively agree with the TB calculations of Reference [23 [281129] , 
i.e., the continuum approach can reproduce all important features of multi-band TB 
calculations. A more quantitative comparison between our results and the TB results 
[271 EHII29] is difficult, partly because the details may depend on the way how the 
material parameters are extracted from the DFT band structure and also because in 
the TB calculations the Zeeman effect was often neglected. 

The LL energies can also be obtained analytically in the approximation where 
and ^ are neglected. We will not show these analytical results here because it turns 
out that an even simpler approximation yields a good agreement with the numerical 
calculations shown in Figures [3 and |2] (see Section 12.2p and offers a suitable starting 
point to develop a theory for the SdH oscillations of the longitudinal conductivity. 

2.2. Approximation of the LLs spectrum 

In zero magnetic held, the trigonal warping term Eq. flHcjl and the third order term 
Eq. fl5^ are important in order to understand the results of recent angle resolved 
photoelectron spectroscopy measurements and in order to obtain a good ht to the DFT 
band structure, respectively [3] . However, as we will show for the calculation of LLs the 
terms Hff and are less important. To see this one can perform another Lowdin- 

partitioning on to obtain effective singe-band Hamiltonians for the VB and the 
CB separately. Keeping only lowest order terms in one hnds that these single-band 
Hamiltonians correspond to a harmonic oscillator Hamiltonian (with different effective 
masses in the VB and CB and for the spin-split bands) and a term which describes a 
linear-in-R 2 splitting of the energies of the LLs in the two valleys. Therefore the LL 
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spectrum can be approximated by 


T^'r,s 

^n,vh 


^vb + 


‘^vb 


TT'US 

^n,ch 


= e 


cb 





+ ^gel-isB^s + r, 

+ T. 


(ба) 

( бб ) 


Here, the following notations are introduced: n = 0,1,2,... is an integer denoting the 
LL index, £yb(cb) ~ ^vb(cb) + 'rAvb(cb)S 2 are the band edge energies in the VB (CB) for a 
given spin-split band s and ~ are cyclotron frequencies. In terms of the 

^ ™-vb’(cb) 

parameters appearing in Eqs. (l2])-(j4l), for r = 1 the effective masses ?7i[b\cb) enter 
the expression of the cyclotron frequencies are given by [3] 




2m 


(1,*) 

vb 




2m 


(1,*) 

cb 





+ /3s + 



(7a) 

{7b) 


where The corresponding expressions for r = — 1 can be easily found 

from the requirement electronic states that are connected by time reversal symmetry 
have the same effective mass. This means that bands corresponding to the same value of 
the product r s have the same effective mass. The third term in Eqs fl6 ap and comes 
from the free-electron term ([2]). The VDB is described by the last term in Eqs. (I6^ . 
(1^ and the valley g-factors are given by 


(s) _ 

yvi,vh ~ f^2 






(8 a) 

( 86 ) 


As one can see from fl8^ - fl^ . depends on the (virtual) inter-band transition matrix 
elements (dg and 7 . Due to the intrinsic spin-orbit coupling, the magnitude of these 
matrix elements is spin-dependent [3]. Note, that g^i is different in the VB and the 
CB. This is in agreement with numerical calculations based on multi-band tight-binding 
models |23129]. For the CB, the details of the derivation that leads to fl^ can be found 
in [3l], for the VB the derivation of (16^ is analogous and therefore it will not be 
detailed here. We note that in variance to Reference [3l], we do not define separately 
an out-of-plane spin g-factor and a spin independent valley g-factor, these two g-factors 
are merged in g)^^ . The response to magnetic held also depends on the free electron 
Zeeman term. The spin-index s to be used in the evaluation of the Zeeman term in 
Eqs. (E3)-(IEQ) follows the spin-polarization of the given spin-split band. For M 0 S 2 , the 
spin polarizations s of each band are shown in Figure |5l other M 0 X 2 (X={S, Se, Te}) 
monolayer TMDCs have the same polarization. For monolayer WX 2 TMDCs the spin 
polarization in the VB is the same as for the M 0 X 2 , but in the CB the polarization of 
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the lower (higher) spin-split band is the opposite [3]. We are mainly interested in how 
the magnetic held breaks the degeneracy of those electronic states which are connected 
by time reversal in the absence of the magnetic held. Using Eqs. the valley 

splitting 5E^b(vb) ~ 9eff ch(vh)f'BBz of these states can be characterized by an ehective 

g-factor 5 'e//,cb(vb) = idesr + where i = 1 ( 2 ) denotes the higher-in-energy 

(lower-in-energy) spin-split band. In the VB the upper index (1) [(2)] is equivalent to 

(t), but in the CB the relation depends on the specihc material being considered 
because the polarisation is diherent for M 0 X 2 and WX 2 materials. Taking hrst the 
M 0 X 2 monolayers one hnds that (see also Figure E]) 

5'e//,vb = {~9e + 9tl,vh) 5'e//,vb = (S'e + 5'Iz,vb) (9®) 

^e//,cb = i9e + 9 I 1 , ch) 9eff, ch = + ^iz,cb) (9^) 

For WX 2 monolayers 5'e/Vvb calculated by (E^, whereas in the CB 

5'e//,cb ~ i~9e + 9tl, ch) ^e//,cb ~ (S'e + 5'Iz,cb) (^9) 

As an example the numerical values of the various g-factors dehned above are given in 
Table [Ufor M 0 S 2 and in Table | 2 ]for WSe 2 . One can see that glfch{vh) ^e comparable 
in magnitude to gg. This explains why the valley splitting is very small for M 0 S 2 in the 
case of the upper spin-split band in the CB (see Figure [I]), whereas the opposite is true 
for WSe 2 (Figure E]). 

As one can see from Fqs. fl 8 ^ and fl^ . g^f depends explicitly on the band-gap 

(s) 

ElJ of a given spin s. In addition, the parameters 7 , as and jdg implicitly also depend on 

is) 

EIJ due to the htting procedure that is used to obtain them from DFT band structure 
calculations [3]. It is known that E^^^ is underestimated in DFT calculations and its exact 
value at the moment is not known for most monolayer TMDCs. Therefore in Reference 
[3] we have obtained two sets of the k • p band structure parameters, the first one using 
EIJ from DFT and the second one using E^J extracted from GW calculations. The 
calculations shown in Figures [Hand [2] were obtained with the former parameter set. As 
shown in Table IH the calculated g-factors depend quite significantly on the choice of the 
parameter set. While there is an uncertainty regarding the magnitude of g)^i ^ we expect 
that the g-factors obtained by using the DFT and the GW parameter sets will bracket 
the actual experimental values. On the other hand, the effective masses are probably 
captured quite well by DFT calculations and therefore the first term in Fqs. fl 6 ^ - fl 6 ^ 
is less affected by the uncertainties of the band structure parameters. The calculations 
in Figures [Hand [2] correspond to the “DFT” parameter set in Tables [H and El 

In order to see the accuracy of the approximation introduced in Fq. in 

Figure [3] we compare the LL spectrum obtained in this approximation and calculated 
numerically using the Hamiltonian ([ 1 ]). As one can see the approximation is very good 
both in the VB and in the GB up to magnetic fields < 20T. For larger magnetic fields 
and large LL indices {n > 7) deviations start to appear between the full quantum results 
and the approximation. The deviations are stronger in the VB which we attribute to 
the larger trigonal warping [3] of the band structure in the VB. To our knowledge the 
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Table 1. Valley g-factors in M 0 S 2 . In the first row the g-factors are obtained with 
the help of DFT band gap, in the second row the g-factors are calculated with a band 
gap taken from the GW calculations. 



Tpi- 

^bg 

< 

9vl,vh 

t 

9vl,vh 

9eff,vh 

(2) 

9eff,.b 

9vl,ch 

t 

9vl,ch 

9eff,ch 

(2) 

9eff,ch 

DFT 

1.66 eV“ 

1.838 eV“ 

0.98 

0.96 

-1.02 

2.96 

- 2.11 

-2.05 

-0.05 

-4.11 

GW 

2.8 eV'' 

2.978 eV^ 

2.57 

2.38 

0.57 

4.38 

-0.52 

-0.6 

1.4 

-2.52 


“ adapted from Reference [3]. 

^ adapted from Reference [55] . 


Table 2. Valley g-factors in WSe 2 . In the first row the g-factors are obtained with 
the help of DFT band gap, in the second row the g-factors are calculated with a band 
gap taken from the GW calculations. 



< 

< 

4- 

9vl,vh 

t 

9vl,vh 

ye//,vb 

(2) 

9eff,vh 

9vl,ch 

t 

^e//,cb 

(2) 

9eff,ch 

DFT 

1.337eV“ 

1.766eV“ 

-0.38 

-0.23 

-2.38 

1.77 

-2.71 

-2.81 

-4.71 

-0.81 

GW 

2.457eV^ 

2.886eV'' 

2.55 

1.9 

0.55 

3.9 

-0.67 

0.13 

-2.67 

2.13 


“ adapted from Reference [3]. 

^ adapted from Reference [36] . 


effects of the non-parabolicity of the band-dispersion on the LL spectrnm has not been 
discussed before for monolayer TMDCs. 



Figure 3. Comparison of the LL spectrum in M 0 S 2 obtained from the two-band model 
and from the single band model, a) The numerically calculated LLs using ([1]) for the 
T = l,s = — lin the VB (squares) and the approximation (solid lines) for LL 
indices n = 0 ... 9. b) the same as in a) but for the for the r = I, s = — I band in the 
CB (squares) and the approximation (l6B (solid lines). 

Given the noticeable uncertainty regarding the exact values of the effective g-factors, 
one may ask which features of the LL spectrum are affected or remain qualitatively the 
same. Looking at Tables [D and E] one can see that in some cases only the magnitude of 
an effective g-factor changes, in other cases both the magnitude and the sign. Firstly, 
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we consider a case which illustrates possible effects of the uncertainty in the magnitude 
of an effective g-factor. In Figure 0] we show the LLs in the lower spin-split CB in M 0 S 2 



Figure 4. Comparison of the LL spectrum in in the lower-in-energy spin-split CB of 
f2l (‘ 2 ') 

M 0 S 2 obtained with a) = —4.11 and b) = —2.52. LLs in different valleys 

are denoted by different colors. 


for the two different 5'e//cb given in Table [T] One can see that in Figure HJ^a) the VDB 
is small, except for the lowest LL, which is clearly separated from the other LLs. If one 
assumes that the LLs acquire a finite broadening then all LLs would appear as doubly 
degenerate except the lowest one in, e.g. an STM measurement. In contrast, the LLs 
are in Figure IH^b) are more evenly spaced and may appear as non-degenerate even if 
they are broadened. 

Secondly, in some cases also the sign of Qef / changes depending on which parameter 
set is used. For g^ff > 0 the LLs in the K valley have higher energy than the LLs in the 
— K valley, while for negative Qeff the opposite is true. We note that in Reference [37] 
Eqs were used to understand the VDB in the excitonic transitions in MoSe 2 . 

The exciton valley g-factor gvi,exc was obtained by considering the energy difference 
between the lowermost LL in the CB and the uppermost LL in the VB in each valley: 


9 ex,vllJ‘BBz — -^n=0,vb) i^n=- 


t=i,4 


r'"'"-“lit p''"——i,t\ 

=0,cb ■^n=0,vb/' 


( 11 ) 


Using Eqs. (I7al) - fl^ . one can easily show that in this approximation the exciton valley 
g-factor is independent of the band gap and it can be expressed in terms of the effective 
masses in the CB and VB [ST] [38] : 


9 ex,vl 4 2 


rrip 


nip 


^cb \^vb\ 


( 12 ) 


Therefore, albeit the effective g-factors in the CB and VB separately are affected by 
uncertainties, the exciton g-factor, in principle, can be calculated more precisely so long 
the effective masses are captured accurately by DFT calculations. The comparison of 
DFT results and ARPES measurements [3| suggest that the DFT effective masses in the 
VB match the experimental results quite well. At the moment, however, it is unclear 
how accurate are the DFT effective masses in the CB. 

Finally, we make the following brief comments. 
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i) In the gapped-graphene approximation, i.e., if one neglects the free electron term 

and the terms ~ in Eqs. fl7^ - fl7^ and in then the lowest LL in the 

CB and the highest one in the VB will be non-degenerate, bnt for all other LLs the 
valley degeneracy wonld not be lifted [H] dne to a cancellation effect between the 
Erst and last terms in Eqs. flb^ and 

ii) By measuring the valley g-factors and the effective masses one can deduce the 
Diracness of the spectrum [18], i.e., the relative importance of the off-diagonal and 
diagonal terms in fl5^ and (j56j), respectively. 

3. Shubnikov-de Haas oscillations of longitudinal conductivity 

As we will show, the results of the Section 12.21 provide a convenient starting point 

for the calculation of the SdH oscillations of the magnetoconductance. Our main 

motivation to consider this problem comes from the recent experimental observation 

of SdH oscillations in monolayer |T8] and few-layer |THl 120] samples. Regarding previous 

theoretical works on magnetotransport in TMDCs, quantum corrections to the low-field 

magneto-conductance were studied in References [25l|26|. A different approach, namely, 

the Adams-Holstein cyclotron-orbit migration theory [39], was used in Reference [ID] 

to calculate the longitudinal magnetoconductance cx^x- This theory is applicable if the 

cyclotron frequency is much larger than the average scattering rate l/fgc- By using 

the effective mass obtained from DFT calculations [3] and taking the measured values 

2 - 

of the zero held electron mobility /ig = ^ and the electron density rig given in 
Reference [H] for monolayer M 0 S 2 , a rough estimate for f^g can be obtained. This 
shows that for magnetic fields R < 15 T the samples are in the limit of oJch^sc ^ 1 and 
therefore the Adams-Holstein approach cannot be used to describe a^x- Therefore we 
will extend the approach of Ando [12] to calculate in monolayer TMDCs because it 
can offer a more direct comparison to existing experimental results. 
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Figure 5. Schematics of the dispersion in the VB and in the CB around the K and 
—K points of the band structure. The spin-split bands are denoted by red and blue 
lines, different colours indicate different spin-polarization. The arrows show the spin- 
polarization for M 0 S 2 . For typical values of doping, the Fermi-level Ep (denoted by a 
dashed line) would intersect only the upper spin-split band in the VB or both spin-split 
bands in the CB. The index (1) and (2) denote the upper and lower spin-split band. 


Before presenting the detailed theory of SdH oscillations we qualitatively discuss 
















LLs and SdH oscillations in monolayer TMDCs 


12 


the role of the doping and the assumptions that we will use. The most likely scenarios 
in the VB and the CB are shown in Figure 0 (a) and (b), respectively. Considering hrst 
the CB, for electron densities ~ 10^^/cm^ measured in Reference [H] both the upper 
and lower spin-split bands would be occupied. In contrast, due to the much larger spin¬ 
splitting, for hole doped samples Ep would typically intersect only the upper spin-split 
VB. Such a situation may also occur for n-doped samples in those monolayer TMDCs 
where the spin-splitting in the CB is much larger than in M 0 S 2 , e.g., in MoTe 2 or WSe 2 . 
For strong doping other extrema in the VB and CB, such as the F and Q points may 
also play a role, this will be briefly discussed at the end of this section. 

We will have two main assumptions in the following. The hrst one is that one 
can neglect inter-valley scattering and also intra-valley scattering between the spin- 
split bands. Clearly, this is a simplihed model whose validity needs to be checked 
against experiments. One can argue that in the VB (see Figure [5](a)) in the absence of 
magnetic impurities the inter-valley scattering should be strongly suppressed because it 
would also require a simultaneous spin-hip. A recent scanning-tunneling experiment in 
monolayer WSe 2 [S] indeed seems to show a strong supression of inter-valley scattering. 
In the CB, for the case shown in Figure |5](b), the inter-valley scattering is not forbidden 
by spin selection rules. Even ii Ep was smaller, such that only one of the spin-split 
bands is populated in a given valley, the inter-valley scattering would not be completely 
suppressed because the bands are broadened by disorder which can be comparable to 
the spin-splitting 2Acb (2Acb = 3 meV for M 0 S 2 and 20 — 30 meV for other monolayer 
TMDCs.) On the other hand, the intra-valley scattering between the spin-split bands in 
the CB should be absent due to the specihc form of the intrinsic SOC, see Eq. (j3]). We 
note that strictly speaking any type of perturbation which breaks the mirror symmetry of 
the lattice, such as a substrate or certain type of point defects (e.g., sulphur vacancies) 
would (locally) lead to a Rashba type SOC and hence induce intra-valley coupling 
between the spin-split bands. It is not known how effective is this mechanism, in the 
present study we neglect it. The second assumption is that we only consider the effect 
of short range scatterers. This assumption is widely used in the interpretation of SdH 
oscillations as it facilitates to obtain analytical results |12]. We note that according to 
References [HlEl], some evidence for the presence of short range scatterers in monolayer 
M 0 S 2 has indeed been recently found. While short-range scatterers can, in general, 
cause inter-valley scattering, on the merit of its simplicity as a minimal model we only 
take into account intra-valley intra-band scattering. 

Using these assumptions it is straightforward to extend the theory of Ando |12] to 
the SdH oscillations of monolayer TMDCs. Namely, as it has been shown in Section 
12.21 for not too large magnetic fields the LLs in a given band can be described by a 
formula which is the same as for a simple parabolic band except that it contains a 
term which describes a linear-in-magnetic held valley-splitting. Then, because of the 
assumption that one can neglect inter-valley and intra-valley inter-band scattering, the 
total conductance will be the sum of the conductances of individual bands with valley 
and spin indices r, s. This simple model allows us to focus on the effects of intrinsic 
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SOC and valley splitting on the SdH oscillations, which is our main interest here. 

Following Reference |12], we treat impurity scattering in the self-consistent Born 
approximation (SCBA) and use the Kubo-formalism to calculate the longitudinal 
conductivity axx (for a recent discussion see, e.g., [IHlIll]). Assuming a random disorder 
potential R(r) with short range correlations (R(r)l/(r')) = Asc5(r — rOi fhe self-energy 
in a given band (r, s) does not depend on the LL index n. It is given 
by the implicit equation 




+ = 


27r/^ 


E 

n=0 


E-E, 


T,S 


(EP* + *E") 


(13) 


where El^^ is given by Eqs (jHaj)- 


The term A^c/27r/^ on the right-hand side of Eq. flT^ 


can be rewritten as 






^ where l/rf} = XsM^/h^ is the scattering rate 
calculated in the Born-approximation in zero magnetic held. As in Section 12.21 the 
upper index i = 1(2) refers to the higher(lower)-in-energy spin-split band in a given 
valley (see also Figure |5]). Using the Kubo-formalism the conductivity coming from a 
single valley and band is calculated as 






dE 


df{E) \ 

I ^xx 


dE 


(E) 


(14) 


where f{E) is the Fermi function and 


(J-r 


■(£) = (MP" E(" + l)R«|Gy(". E)G"{n +!,£)- GY{n, E)GJ(n + 1, B)](15) 


n=0 


Here G^^{n,E) and G^J^{n,E) are the retarded and advanced Greens-functions, 
respectively. Vertex corrections are neglected in this approximation. Since we 
neglect inter-valley and intra-valley inter-band scattering, the disorder-averaged Greens- 
function E) = [E — Ei^’^ — is diagonal in the indices r, s and in the LL 

representation it is also diagonal in the LL index n. The total conductivity is then 
given by axx = s '^xx where the summation runs over occupied subbands for a given 
total electron (hole) density Ue (uh). In general, one has to determine -|- ^E[’* by 
soving Eq. flT^ numerically. The Greens-functions can be then calculated and a^’^ 
follows from Eq. (lT5ll . It can be seen from Eq. flTdll that at zero temperature E'^’^(E) 
and a'^’^{E) has to be evaluated at E = Ep. In the semiclassical limit, when there are 
many occupied LLs below Ep, i.e., Hooc'^ Ep, one can derive an analytical result for 
a^’^, see References [l2l HS] for the details of this calculation. Here we only give the 
hnal form of axx and compare it to the results of numerical calculations. 

As mentioned above, the situation depicted in Figure |5](a), i.e., when there is only 
one occupied subband in each of the valleys is probably most relevant for p-doped 
samples. One hnds that in this case the longitudinal conductance is 


(^xxj n’o — 


1 + (^•lbT-ic^)2 


1 + (^lbT-iP)2 
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9effiapB^ 4riP)^ 
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2 ( 1 ) 2 ( 1 ) 

Here (Tq = Ep = is the zero field conductivity per single valley and band, 

^vb 

Uh is the total charge density and we assumed Ep. The amplitudes Ai ^2 

and B are given by 


TT 


= COS I 


m. 


( 1 )' 

vb 




B = 


27i^kBT/hw^^ 
sinh (27r‘^kBT/huj^') 


A 2 = sin 


(^ ( 1 ) 

1 2 ^e//,vb 



(17a) 

(17(,) 


where is the Boltzmann constant and T is the temperature. One can see that 
Eqs. flT6l) -fll7fe|l are very similar to the well known expression derived by Ando m 
for a two-dimensional electron gas (2DEG). The valley-splitting, which leads to the 
appearance of the amplitudes ^ 1 , 2 , plays an analogous role to the Zeeman spin-splitting 
in 2DEG. Therefore, under the assumption we made above, the uncertainty regarding 
the value of the effective g-factors affects the amplitude of the oscillations but not their 
phase. The term proportional to fiBBz/Ep in Eq. fllbp is usually much smaller than the 
first term. Thus, it can be neglected in the calculation of the total conductance, but 
may be important if one is interested only in the oscillatory part of axx, see below. 

We emphasize that Eq. ffT6|) is only accurate if <C Ep. However, in 

semiconductors, especially at relatively low doping, one can reach magnetic field values 
where the cyclotron energy becomes comparable to Ep. In this case the numerically 
calculated axx may differ from Eq. flT^ 0. It is known that, e.g., WSe 2 can be relatively 
easily gated into the VB, and a decent Hall mobility was recently demonstrated in few- 
layer samples in Reference [15]. As a concrete example we take the following values [T5] : 
Uh = —4 * 10^^/cm^ and Hall mobility hh = 700cm^/Vs. By taking = —0.36me 
[3] the Fermi energy is Ep ^ —26.6meV and using that fip = we obtain 

ric^ = 1.4 X 10“^^s. The amplitude of the oscillations should become discernible when 
= finBz < 1, i.e., for magnetic fields Bz > lOT, while at Bz = 14.28T, which 
corresponds to ~ 1, there are around six occupied LLs. One can expect that for 

Rz < 15T the LL spectrum is well described by Eq fl6^ . however, since the number of 
LLs is relatively low, there might be deviations between the analytically and numerically 
calculated axx- In Figure |6]^a) we show a comparison between the analytical result 
Eq. (1T6|) and the numerically calculated longitudinal conductance at zero temperature. 
The effective g-factors used in these calculations are given in Table |2l 

One can see that for larger magnetic fields the amplitude of the oscillations is 
not captured very precisely by Eq. ffT6|) but the overall agreement with the numerical 
results is good. Next, in Figures |6]^b) and (c) we compare the oscillatory parts axx,osc 
of the longitudinal conductivity obtained from numerical calculations and from Eq (fT6|) 
using two different fi'e//,vb values. In the case of the numerical calculations axx,osc was 
obtained by subtracting the smooth function 2/(1 -|- from axx- According to 


f From a theoretical point of view, in strong magnetic fields one should also calculate vertex correlations 
to axx, but this is not considered here. 
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Figure 6. Comparison of the numerically (symbols) and analytically (solid line) 
calculated zero temperature conductivity for WSe 2 for the situation depicted in Figure 
[5] a) total conductivity, <7e//vb = 0.55. b) and c) comparison of the oscillatory parts 

of Oxx- In b) 5e//vb = 0.55, while in c) ffe//vb ~ —2.38 (see Tabled. The figures 
correspond to a magnetic field range of about 5.7 — 15.7 T. 


Eq. flT6|) . the valley-splitting of the LLs and the different effective g-factors should only 

affect the amplitude of the oscillations. While the amplitude of the oscillations indeed 

depends on agreement is better for f/ez/vb ~ than for 

5'e//vb ~ —2.38. For the latter case the position of the conductance minimuma start 

to differ for large magnetic held, whereas the maximuma in a^x agree in both hgures. 

These calculations illustrate that Eq. flTHll may not agree with the numerical results 

when there are only a few LLs below Ep. 

We now turn to the case shown in Figure [5](b) when both spin-split subbands are 

populated. The total conductance is given by the sum of the conductances coming 

from the two spin-split subbands : axx = crix + (jfx- Since the effective masses in the 

spin-split bands are, in general, different, the associated scattering times ric^ and Ts^ 

calculated in the Born-approximation are also different. We dehne fgc = 'ric^ + Ts^ and 
2 - 

do = i^^Ep, and obtain for the magnetoconductance 
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In Eq. ffTSj) we have neglected terms which are ~ fiuBz/Ep. The result shown in Eq. fflSj) 
is similar to the multiple subband occupation problem in 2DEG 113 SHI HZ] . The valley 
splitting affects the amplitude of the oscillations (see Eq. (1196^ 1. whereas the intrinsic 
SOC can influence the amplitude of the oscillations [see Eq. (IIQqI ] and it also leads 
to a phase difference [Eq. flTSl) ] between the oscillations coming from the two spin-split 
subbands. 

The situation depicted in Figure[H][b) is easily attained, e.g., in the CB of monolayer 
M 0 S 2 , where DFT calculations predict that the spin-splitting is 2Acb = 3meV and 
therefore both subbands can be populated for relatively low densities. Our choice of 
the parameters for the numerical calculations shown below is motivated by the recent 
experiment of Cui et al. [18], where SdH oscillations in mono- and few layer M 0 S 2 
samples have been measured. We use rig = 10^^/cm^ and mobility fin = lOOOcm^/Vs. 
The effective masses are chosen as = 0.46me, = O.dSmg and the spin-splitting 

in the CB is 2Act = 3meV [3]- Using these parameters we hnd Ep = 28.43 meV. Since 
the effective masses are rather similar, the scattering times calculated from are close 
to each other: ric^ ~ ric^ 2.6 x 10“^^s, i.e., they are almost twice as long as in the 
case of WSe 2 . The oscillations in axx should become discernible for > 7T, and at 
Bz = lOT there are ten LLs in both the lower and the upper spin-split CB in each 
valley. We will focus on the oscillatory part axx,osc = o'is,osc + o'^x,osc of the conductance, 
since this contains information about the spin and valley splittings. As in the previous 
example, we hrst calculate ctxx numerically using Eqs. (1 13 p - (1 15 p and obtain (Txx ,osc by 
subtracting the smooth function 2C^^y\i -f [ujc'^TscY]. We than compare these results 
to the oscillations that can be obtained from Eq. flTSl) . 

In Figures [7](a) and (b) we show the numerically calculated and aix for the 
two different sets of g-factor values given in Table [1] as a function of cJcAc which 
was introduced as a dimensionless scale of the magnetic held. Here Ac = with 

^ = \lfsc = \/Tsl^Tsc^■ All calculations are at zero temperature. One can 
observe that due to the CB spin splitting 2Act the oscillations of and will not be 
in-phase for larger magnetic held. This ehect is expected to be even more important for 
TMDCs having larger 2Acb than M 0 S 2 and leads to more complex oscillatory features 
in the total conductance than in the previous example of p-doped WSe 2 where 
only one band in each valley contributed to the conductance. One can also observe 
that in Figure [7j[b) additional peaks with smaller amplitude appear in a^J^osc for larger 
magnetic helds, while there are no such peaks in a^J^osc in Figure [7j[a). The origin of 
this behaviour can be traced back to the diherent valley-splitting patterns shown in 
Figure m The valley splitting of the LLs in Figure llj[a) is small (except for the lowest 
LL), while in Figure |l][b) all LLs belonging to diherent valleys are well-separated for 
larger helds and this leads to the appearance of the additional, smaller amplitude peaks 
in a^J^osc in Figure [7j[b). The comparison between the numerically calculated total 
oscillatory part + o'^x,osc and the corresponding analytical result given in 

Eq. flTSll is shown in Figures [7j[c) and (d). The agreement between the two approaches 
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Figure 7. Oscillations of Uxx in n-doped M 0 S 2 . a) numerically calculated aiH 
(red squares) and aisl (blue circles) using ~ —0.05 and 5e//cb = —4.11. b) 

numerically calculated aill (red squares) and CTosc (blue circles) using Sg^j^b — 

and g^ffch — —2.55. c) The total oscillatory conductance CTosc = Cosl + CTosc 
corresponding to a) (red squares) and the analytical result calculated from Eq. (fT^ 
(solid line), d) the same as in c) but corresponding to b).The figures correspond to a 
magnetic held range of about 4 — 14 T. 


is qualitatively good for < 1. However, for larger magnetic fields where uicfg^ > 1 
the amplitude of the oscillations start to differ. In this regime the oscillatory behaviour 
in can be quite complex, influenced by both the valley splitting and also by the 

intrinsic SOC splitting of the bands. 

We have tried to analyze the experimental results by Cui et al. [H] using the 
theoretical approach outlined above. To this end we have hrst calculated crxx,exp{Bz) by 
inverting the experimentally obtained resistance matrix and normalized it by the zero- 
field conductance o'xx,exp{0)- To simplify the ananlysis, we assumed that the effective 
masses are the same in the two spin-split CB: = OASnie, and hence ric^ = 

ric^ We then htted (Txx,exp{Bz)/(T xx{0) by the function fo{Bz) = C + A/{1 + {yiqBzf), 
where the amplitudes A, C and the quantum mobility Hq are £t parameters. This 
function, according to Eq. (lT8l) . should give the smooth part of the conductance. The 
fit was performed in the magnetic field range [4T — 15T]: for smaller fields the weak- 
localization corrections might be important which are not considered in this work, 
while in larger magnetic held the semiclassical approximation may not be accurate. 
We have found that a^x^exp can be approximated quite well by /o(H^)- The most 
important parameter that can be extracted from the £t is the quantum scattering 
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time Tsc,q-, which is obtained from Tsc,q = We hnd that it is roughly 3.5 times 

shorter than the transport scattering time Tsc,tr that follows from the measured Hall 
mobility i^h = lOOOcm^/Vs. The ratio Tsc,tr/Tsc,q depends to some extend on the 
fitting range that is used, but typically it is Tsc,tr/Tsc,q > 2. This difference may be 
explained by the fact that small-angle scattering is unimportant for Tsc,tr but it can 
affect Tsc,q- We note that Cui et. al. [18] has also found that the Tsc,tr is larger 
than Tsc,q, but they have used the amplitude of the longitudinal resistance oscillations 
in the magnetic field range 10 — 25T to extract and obtained Tsc,trl'^sc,q ~ 1-5. 
The significantly shorter Tgc.q makes it difficult to analyze the magnetic oscillations 
in a quantitative way using Eq. ffTSj) . Namely, it implies that oscillations should be 
discernible for > 15T, i.e., for magnetic fields where only a few LLs are occupied and 
the semiclassical approximation may not be accurate. Using fo{Bz) we then extracted 



Figure 8. Comparison of the theoretical results and the measurements of Cui et 
al. [18]. The measured conductance oscillations cr^x^xp (squares) and the fitting of 
the function fosdBz) (see Eq. (l20l) l using Ue = 1.1 * 10^^/cm^, 2Acb = 3meV (blue) 
He = 1.31 * 10^^/cm^, 2Acb = 3meV (black), and Ug = 1.31 * 10^^/cm^, 2Acb = 5 meV 
(purple). 


the oscillatory part o'xx,osc{d3z) = o'xx,exp{Bz) — /o(-B^) of the conductance and fitted it 
with the function 


fosc{Bz) — —A 
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( 20 ) 


where Di ^2 are fitting parameters. As one can see in Figure [HI the fit can qualitatively 
reproduce the meauserements, but the complex oscillations between 15 — 22T are not 
captured. We also note that a somewhat better fit can be obtained if we assume that 
the charge density is larger than what is deduced from the classical Hall measurements 
(see the black line in Figure IH]) and if we choose the spin-splitting larger than the value 
obtained from DFT calculations (purple line). In all cases we find, however, that the 
fit parameters Di and D 2 differ quite significantly in their magnitude, which is difficult 
to interpret in the present theoretical framework. This might indicate that additional 
scattering channels, such as inter-valley scattering, would have to be taken into account 
for a quantitative theory. 
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Finally, we would briefly comment on the relevance of the other valleys in the band 
structure for the SdH oscillations. Regarding p-doped samples, the F point might, in 
principle, be important for M 0 S 2 . However, according to DFT calculations [3] and 
ARPES measurements [49] the effective mass at the F point is significantly larger than 
at duK and therefore we do not expect that states at F would lead to additional SdH 
oscillations. Nevertheless, they can be important for the level broadening of the sates 
at duK because scattering from ±iF to F does not require a spin-flip [3]. In the case of 
other monolayer TMDCs the F valley is most likely too far away in energy from the top 
of the VB at ±iF to influence the transport for realistic dopings [3|. The situation can 
be more complicated for n-doped samples, especially for WS 2 and WSe 2 . For these two 
materials the states in the six Q valleys are likely to be nearly degenerate with the states 
in the dcK valleys. Therefore the Q valleys might be relatively easily populated for finite 
n-doping and, in contrast to the F point, the effective mass is comparable to that in the 
iA' [3] valleys. Therefore they may contribute to the SdH oscillations. They would also 
affect the level broadening of the dcK valley states because scattering from K {—K) to 
three of the six Q valleys is not forbidden by spin selection rules [3]. Furthermore, we 
note that in the absence of a magnetic held the six Q valleys are pairwise connected 
by time reversal symmetry. Therefore, taking into account only the lowest-in-energy 
spin-split band in the Q valleys, the LLs belonging to the Q valleys will be three-fold 
degenerate: the magnetic held, similarly to the case of the K and —K points, would lift 
the six-fold valley-degeneracy. The effective valley g-factors, however, might be rather 
different from the ones in the ±iF valleys. For n-doped monolayer M 0 X 2 materials the 
situation is probably less complicated because the Q valleys are higher in energy and are 
not as easily populated as for the WX 2 monolayers. For M 0 S 2 monolayers, therefore, 
one can neglect the Q valleys in first approximation. 

4. Summary 

In summary, we have studied the LL spectrum of monolayer TMDCs in a k • p theory 
framework. We have shown that in a wide magnetic held range the effects of the 
trigonal warping in the band structure are not very important for the LL spectrum. 
Therefore the LL spectrum can be approximated by a harmonic oscillator spectrum 
and a linear-in-magntic held term which describes the VDB. This approximation and 
the assumption that only intra-valley intra-band scattering is relevant allowed us to 
extend previous theoretical work on SdH oscillations to the case of monolayer TMDCs. 
In the semiclassical limit, where analytical calculations are possible, it is found tha 
the VDB affects the amplitude of the SdH oscillations, whereas the spin-splitting of 
the bands leads to a phase difference in the oscillatory components. Since in actual 
experimental situations there might be only a few occupied LL below Ep, we have 
also performed numerical calculations for the conductance oscillations and compared 
them to the analytical results. As it can be expected, if there are only a few LLs 
populated the amplitude of the SdH oscillations obtained in the semiclassical limit does 
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not agree very well with the results of numerical calculations. This should be taken 
into account in the analysis of the experimental measurements. We used our theoretical 
results to analyze the measured SdH oscillations of Reference [H] . It is found that the 
quantum scattering time relevant for the SdH oscillations is signihcantly shorter than 
the transport scattering time that can be extracted from the Hall mobility. Finally, we 
briefly discussed the effect of other valleys in the band structure on the SdH oscillations. 
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